(* ::Package:: *)

SetDirectory["C:\\Projects\\branches\\HMC_core\\output"]


param=OpenRead["..\\input\\params_hmc.txt"];
Skip[param,{Number}];
Tstart=Read[param,Number];
Tstep=Read[param,Number];
Tend=Read[param,Number];
Pacc=Read[param,Number];
Nmd=Read[param,Number];
Nmc=Read[param,Number];
Skip[param,{Number}];
Nruns=Read[param,Number];
NumDim=Read[param,Number];

Close[param];

templist=Range[Tstart,Tend,Tstep];

funcname="Griewank\\Mod";


sdOfEpotEnd=Array[c,{Length[templist],Nruns,NumDim}];
sdFinal=Array[b,{Length[templist],2}];
rejected=Array[rej,{Nruns*Length[templist]}];
plot=Array[plott,{NumDim}];
xsum=Array[0&,NumDim];
xsqsum=Array[0&,NumDim];
y=Array[0&,NumDim];

(*whichVar=1;*)


(*Initial Config*)

ifrej=OpenRead["nRejected.reject"];

For[i=1,i<=(Nruns*Length[templist]),i++,
 Skip[ifrej,{Number}];
 rejected[[i]]=Read[ifrej,Number];
];

Close[ifrej];


For[whichVar=1,whichVar<=NumDim,whichVar++,
 Print["whichVar = ",whichVar];
 index=1;
 For[temp=Tstart,temp<=Tend,temp+=Tstep,
  (*Print["Index = ",index, "\nTemp = ", temp];*)
  For[iNruns=1,iNruns<=Nruns,iNruns++,
   whichFile=(Nruns*(index-1)+iNruns);
   loadstring=".\\"<>funcname<>"\\HMC_coord"<>ToString[whichFile]<>".coord";

   fid=OpenRead[loadstring];
   numPoints=0;

   For[iNmc=1,iNmc<=((Nmc+1)-rejected[[whichFile]]),iNmc++,
    Skip[fid,{Number}];
    For[iNumDim=1,iNumDim<=NumDim,iNumDim++,

     x=Read[fid,Number];
     xsum[[iNumDim]]+=x;
     xsqsum[[iNumDim]]+=x^2;

    ];

   ];

   Close[fid];

   For[iCoord=1,iCoord<=NumDim,iCoord++,
    sdOfEpotEnd[[index,iNruns,iCoord]]=Sqrt[(xsqsum[[iCoord]]/((Nmc+1)-rejected[[whichFile]]))-((xsum[[iCoord]]/((Nmc+1)-rejected[[whichFile]]))^2)];
   ];
   Clear[xsum,xsqsum,x,y,numPoints];
   xsum=Array[0&,NumDim];
   xsqsum=Array[0&,NumDim];
   y=Array[0&,NumDim];

  ];

  sdFinal[[index,2]]=Sum[sdOfEpotEnd[[index,j,whichVar]],{j,1,Nruns}]/Nruns;
  sdFinal[[index,1]]=templist[[index]];
  index+=1;
 ];

 plot[[whichVar]]=ListPlot[{Tooltip[Re[sdFinal]],Table[{m,(5*Sqrt[3])/2},{m,Tstart,Tend,Tstep}]},PlotRange->{Full,Full},Joined->{False,True,True},AxesLabel->{Temperature,SD},AxesOrigin->{Tstart,0},PlotStyle->PointSize[Medium]];
 Clear[sdFinal];
 sdFinal=Array[b,{Length[templist],2}];
];


{plot[[1]],plot[[NumDim]]}


Gn=1+(((Subscript[x,1]^2)+(Subscript[x,2]^2))/4000)-Cos[Subscript[x,1]]Cos[Subscript[x,2]/Sqrt[2]]


str=OpenRead["C:\\Projects\\branches\\HMC_core\\output\\Griewank\\Mod\\HMC_coord30.coord"];

a=Array[aa,{((Nmc+1)-rejected[[30]]),2}];
Skip[str,{Number}];
For[i=1,i<=((Nmc+1)-rejected[[30]]),i++,
 a[[i]]=Read[str,{Number,Number}];
 Skip[str,{Number}];
];

ContourPlot[Gn,{Subscript[x,1],-100,100},{Subscript[x,2],-100,100},AxesLabel->Automatic,ContourShading->None,Epilog->{Black,Line[a],Black,Point[a]}]


Plot3D[Gn,{Subscript[x,1],-87,-57},{Subscript[x,2],14,16},AxesLabel->Automatic,PlotRange->Full]


N[5Pi]


Subscript[x,2]=2.2;


Unset[Subscript[x,2]]


Plot[Gn,{Subscript[x,1],-87,-57},AxesLabel->Automatic,AxesOrigin->{-40,0},PlotRange->{0,10}]


tGn=-Cos[x]Cos[y/50]


N[Cos[67/50]]


Plot3D[tGn,{x,-87,-57},{y,-100,100}]


y=0


Plot[tGn,{x,-100,100},PlotRange->{-1,1}]


Unset[y]


Close[str]
